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ABSTRACT 

We present a new model for the spectral evolution of Pulsar Wind Nebulae inside 
Supernova Remnants. The model couples the long-term dynamics of these systems, as 
derived in the 1-D approximation, with a 1-zone description of the spectral evolution 
of the emitting plasma. Our goal is to provide a simplified theoretical description 
that can be used as a tool to put constraints on unknown properties of PWN-SNR 
systems: a piece of work that is preliminary to any more accurate and sophisticated 
modeling. In the present paper we apply the newly developed model to a few objects 
of different ages and luminosities. We find that an injection spectrum in the form of 
a broken-power law gives a satisfactory description of the emission for all the systems 
we consider. More surprisingly, we also find that the intrinsic spectral break turns out 
to be at a similar energy for all sources, in spite of the differences mentioned above. 
We discuss the implications of our findings on the workings of pulsar magnetospheres, 
pair multiplicity and on the particle acceleration mechanism(s) that might be at work 
at the pulsar wind termination shock. 

Key words: supernova remnants; pulsars: general; stars: winds, outfiows; radiation 
mechanisms: non-thermal 
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1 INTRODUCTION 

Pulsars (PSRs) are known to release most of their ro- 
tational energy in the form of a relativistic magnetized 
wind, whose particle component is mostly made of electron- 
positron pairs, with, possibly, a minority of ions, while the 
magnetic field is a lmost purely toroid a l far away from the 
light cylinder (e.g. iGoldreich fc Julian! (|l969l )). The wind is 
initially confined by the slowly expanding Supernova Rem- 
nant (SNR), which leads to the formation of a termination 
shock at some distance from the pulsar, and of a bubble of 
relativistically hot fluid beyond it. This is the so-called Pul- 
sar Wind Nebula (PWN), that shines through synchrotron 
and Inverse- Compton emission from radio wavelengths up 
to 7-rays. 

The typical energy that a pulsar injects into the PWN 
during its entire lifetime is about 10~^ of the typical energy 
of the SN explosion (10""^ erg). Therefore, the presence of 
an energetic PSR, has little effect on the global evolution 
of the SNR. In contrast, the evolution of the PWN strongly 
depends on its interaction with the SNR (see Sec. [3]). 
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Modeling of the evolution and dynamical properties of 
PWN-SNR systems has much improve d since the pio n eering 
one -dimensional analyt i cal wo rks by iRees fc GunrJ (jlQTJ) 
and Kennel fc Coronitil (|l984a). In recent years, the impres- 
sive advances in X-ray observations, due to CHANDRA and 
XMM-Newton, have led to a renewed interest in these ob- 
jects. A major effort has been devoted to explaining the ob- 
served X-ray properties of the Crab Nebula and other young 
systems through multi-di mensional time-dependent ap- 
proaches to their dynamics IJKomissarov fc Lyubarskvll2004l: 



Del Zanna et al.l2004l.l2006l) Several authors JBlondin et all 



20011: 



der Swal uw et al.|[200ll : iBucciantini et al.ll2004bl : 
de Jager et al. 200g) have presented numerical evolution- 
ary models, extended to advanced ages (typically up to 
• ^ 10^ yr), and to the case of fast moving pulsars 
(|van der Swaluw et al.ll2004l '). but in fact, fully multidimen- 
sional models are still lacking in these cases. The reason for 
this is that detailed modeling of older systems is far more 
difficult than for the young ones: the latter are brighter, 
hence providing us with high quality multi-frequency data, 
and are often better constrained in terms of the most rel- 
evant physical quantities, such as age, distance and pulsar 
energetics. 
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The overall evolution of old objects is much more com- 
plex than the smooth expansion appropriate to describe 
young systems. In general, even within the most simplistic 
approach, the long-term evolution of a PWN-SNR system 
depends on many parameters: the SN energy, the mass of 
the ejecta, their density structure, the d ensity of the ISM, 
the PS R luminosity, its spin-down history (jBucciantini et alj 
[2004bj). A model for particle injection and spectral evolution 
must be added to all of this, if the aim is that of deriv- 
ing information from direct comparison with observations. 
At present, multi-dimensional numerical simulations are not 
suitable to properly investigate the large parameter space 
that older systems may span. 

The situation is actually delicate also for many young 
systems. Indeed, accurate multidimensional results, includ- 
ing emission maps to be directly compared with observa- 
tions, have been so far presented only for the Crab Neb- 
ula. But this is a very special object, for which measure- 
ments are available through almost the entire electromag- 
netic spectrum and almost all the aforementioned parame- 
ters are known, including the age and the pulsar rotational 
history (initial spin frequency and braking index). For other 
young systems, like 3C58 or MSH 15-52, the quantity and 
quality of information available is much poorer than in the 
case of Crab and simplified evolutionary models are again 
the only possible tool for investigation, at least for a start. 

Models for the evolution of the emission proper- 
ties of PWNe h ave a long history, from the work by 
[Pacini & Salvati fl973), to t he most recent developments 
byjdc Jagcr ct al. (200i), and ICelfand et all ([2009). These 
later works have attempted at taking into account both 
the dynamics of PWN-SNR systems, as derived from ac- 
curate numerical simulations, and the spectral evolution of 
the particle di stribution function. In particular, the work by 
ICelfand et al., (2009) is exemplary for showing the richness 
of behaviors that is found even within the one-zone treat- 
ment, and the degree of complexity that one can investigate 
using simple tools. 

While constraining the unknown parameters is prelimi- 
nary to detailed multi-dimensional modeling of the systems, 
which can only sample a much more limited range of parame- 
ters, the question arises of how strong and reliable the results 
obtained with a 1-zone approach are. The recent success of 
2-d models of PWNe at reproducing so many of the ob- 
served features has interesting consequences in this respect, 
in that it lends support to the 1-zone description over the 
1-dimensional approximati on. The latter inevitably leads to 
onion-like structures (Ken nel fc Coronitill 198 4b 1 , where par- 
ticles injected at earlier times are located at larger radii, 
with a one-to-one correspondence between the properties of 
the particle distribution function and the distance from the 
central PSR. Comparison of multidimensional models with 
observations has shown that such ordered structure is not 
what is realized in PWNe, where an overall turbulent flow 
is injected at the terminat ion shock and sur vives through 
most of the nebular extent ([Camus et al.ll2009l ). For the typ- 
ical flow times, that are much smaller than the PWN age, 
this turbulent flow-structure leads to efficient mixing, which 
is what a 1-zone description is based on. This also finds 
some direct hints from observations where detailed data are 
available, as in the case of the Crab Nebula: the remark- 
able homogeneity of radio spectral index across the nebula. 



and the extent of the X-ray emitting region which exceeds 
the predictions of any simple 1-dimensional flow model in 
the absence of particle diffusio n far more efficient than the 
Bohm rate (lAmato et al.ll2000l '). 

In the following we present a s implified evolutionar y 
model, similar in spirit to the work bv lCelfand et al.l ([2009|). 
Our approach combines analytical results and numerical 
simulations for the description of the dynamics of the PWN- 
SNR system and the evolution of the particle distribution 
function within the PWN: the key features of the interaction 
between a PWN and the surrounding SNR are adequately 
reproduced during different evolutionary stages. The model 
is here applied to a few systems of various ages, both young 
and older ones with the aim of showing how this kind of 
modeling can help the extraction of information on the sys- 
tem properties from multi-wavelength observations: we try 
to clarify how these properties depend on the different pa- 
rameters, also discussing potential degeneracies and the ob- 
servations that would best serve to disentangle them. 

Special attention is devoted to the particularly interest- 
ing and difficult task of constraining the particle injection 
in PWNe. The relevance of this issue in terms of pulsar 
physics and shock acceleration physics deserves a somewhat 
extended discussion, which will be the topic of next section 
(Sec.[2j. 

The rest of the paper is organized as follows. Sec. |3l 
Sec. 131 and Sec. [S] describe, respectively, how the nebular 
evolution is implemented in our model, the scheme we adopt 
to follow the evolution of the particle distribution function 
and the emission model. In Sec. [S] we show the application 
of our model to several objects, and, finally, in Sec. [7] we 
discuss the implications of our findings. 



2 PARTICLE INJECTION IN PWNE 

Forty years of research on Rotation Powered Pulsars have led 
to a conundrum: observations of young Pulsar Wind Nebu- 
lae (PWNe) have clearly established large particle injection 
rates into the nebulae, well in excess of the electr odynamic 
minimum suggested bv lColdreich fc JulianI (|l969f) . The lat- 
ter corresponds to: 



Nn 






7.6 X 10^ 



Pj'a 4 X 10-13 



1/2 



(1) 



where <I> — \/ Ljc is the magnetospheric potential, L = IQ£l 
is the spin down luminosity, with / the neutron star's mo- 
ment of inertia (745 = 7/10"*^ g cm^), Q, = I-kJP the 
stellar angular velocity and P the rot ation period (P33 = 
P/33 msec). The simple st estimates ([Rees fc GunnI Il974l : 
iKennel fc Coronitilll984al'l s how that X-ray emitting PWNe 
([Kargaltsev fc Pavlovll2008l ) with synchrotron cooling times 
of the X-ray emitting particles well less than the nebular 
ages have electron (and positron) injection rates much larger 
than shown in Eq. ([!}: values reported in the literature for 
the multiplicity k = N /Nr are typically ^^ 10"*, based on 
analysis of the X-ray emission, a result consistent with the- 
oretically derived pair creation rates for the young, high 
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voltage pulsars which have been s u bjected to such analy - 
sis [e.g. iKennel fc Coronitil j 198431 ): iGaensler et~all l|2002l ): 
iHibschman fc AronsI I 2001f l]. Indeed, the fact that the in- 
ferred injection rates are this large is one of the major pillars 
of support for the theoretical conclusion that pulsars have 
substantial outflows of e , this being the only known means 
through which such cool objects can have winds denser than 
the electrodynamic minimum. 

Such high density flows with k 2> 1 support the 
use of the force-free limit of MHD theory in mod- 
eling the torques on rotation pow e red neutron star s 
iKalapotharakos fc Contopoulod l2009l : ISpitkovskvl l2006h . 
As already mentioned, MHD theory has been quite suc- 
cessful in modeling the multidi mensional dynamics and 
appe a ,rance of young PWNe (IKomissarov fc Lyubarskvl 



appe a ,rance ot young rw i Ne (IKomissarov fc Lvu 
12004 iDel Zanna erani2004 l2006l : IVolpi et alj|2008l 



a). Such 



models confirm the earl y inference I Rees fc GunrJ Il974l : 
iKennel fc Coronitilll984al ) that the system behaves as if just 
upstream of the pulsar wind termination shock (TS) the 
plasma has low magnetization a^j = (_B^/87r7,i,n±m±c^) ^ 
1 {n± is the number density of pairs, the particle density is 
2n±): typically, the average of cr,„ in latitude with respect 
to the rotation axis i^ ~ 0.02. 

The nebular dynamics in the MHD model is insensi- 
tive to the specific value of the upstream 4-velocity u^ = 
cj3wjw, as long as 7^, ^ 1. By applying their ID MHD 
model to the optical and harde r photon emission in the 
Crab. IKennel fc Coronitil l|l984bl ) in ferred 7^, ~ 10*^ ■'^ with 
a-w ~ 0.005. Mo re modern models (|Del Zanna et al.l 120061 : 
IVolpi et al.ll2008l ) require latitude averaged CTu, to be some- 
what larger (so that magnetic hoop stress can create the jet 
component of the torus-jet structure), while account of the 
high energy synchrotron emission requires particle spectra 
at the TS with parameters similar to the ID Kennel and 
Coroniti model: again j^ ~ 10®, although this is not explic- 
itly stated since the distribution functions were not tied to 
the specifics of the MHD flow. 

If the apparent low value of cr,„ means that the wind 
just upstream of the TS really is weakly magnetized, a^ 
and 7u, are closely tied to the pair multiplicity k± — k/2, 
since when a <^ 1 the wind carries the rotational energy lost 
from the pulsar in kinetic energy of the flow, L = Mc^jw 
Since M = 2K±m±N{i = 2K±m±c$/e, and L = c<I>^, 
then 7„ = e^/2m±c^K±. In the case of the Crab Neb- 
ula and pulsar, <I> = 4 x 10^® V, while extant theoretical 
mode ls of pair creation, from polar caps, slot gaps or outer 
gaps teibsch man fc Aronsll200ll : IChenJ2007l : lHirotanill2008l : 
[Haxding 2004) all yield k± ~ 10*, for young pulsars with 
$ > lO^^V. Thus theory also says j^ ~ 10® - if one confines 
one's analysis to the high frequency emission from PWNe, 
theory and observation appear to be in good accorqj. 

Using the X-ray emitting particles in PWNe, whose syn- 
chrotron cooling times are short, takes advantage of such 
PWNe being calorimeters for the contemporary particle in- 
jection rates, which gives the results some independence 



^ This analysis has been done in full only for the Crab Nebula, 

however. 

^ We neglec t a possible c omponent of heavy ions in the wind 's 

composition lJGallant fc A rons 1994; Spi tkovskv fc Arong 120041 : 

this subject will be discussed later on 



from the uncertainties of evolutionary models. However, ra- 
dio and infrared emitting particles are much more numer- 
ous than the X-ray emitting ones, due to the rapid de- 
cline of synchrotron emissivity with declining part icle en- 
ergy JKennel fc Coronitil Il984bl : ICallant et fflZ.ll2002l ). Mod- 
els of the underlying pulsar must be able to account for 
all the radiating relativistic particles found in the PWNe 
- the spectral continuity in several systems demonstrated 
in the results described below suggests the lower energy 
particles are indeed injected by the pulsar, rather than 
being accelerated out of the non-relativistic material of- 
ten found embedded within PWNe. Thus measuring the 
full relativistic particle content in PWNe is an impor- 
tant experimental input into modeling of pulsars, distinct 
from the modeling of the energy inpu t into the nebu- 
lae. For the latter the fo rce free model (|Spitkovskvl I2OO6I : 
iKalapotharakos fc Conto poulos 20091) provides the essential 
description, but unfortunately this is independent of the 
multiplicity, so long as k ^ 1. 

We study several young PWNe (Crab, 3C58, B1509, 
Kes75) for which the data can be used to reasonably con- 
strain the pair injection rate - in principle one would like to 
use nebulae with ages known, and reasonably complete (in- 
cluding near and far infrared and millimeter) spectral energy 
distributions (SEDs) are required. We show that all these 
systems have spectral continuity from the radio through the 
infrared to the X-ray bands, suggesting a single source for 
the radiating particles is present in each system. We use 
these full SEDs to derive new estimates for the pair cre- 
ation multiplicity. We discuss the constraints these results 
set on PSR pair production gap models - Polar Caps, Slot 
Gaps, and Outer Gaps - concluding that no existing model 
adequately explains particle injection rates. 

We also discuss alternate hypotheses, that low energy 
particles are picked up from thermal gas in the nebulae, or 
are fossils left over from some unna med acceleration pro- 
cess i n the early history of the nebula JAtovan fc Aharoniani 
Il996l ) , or represent the effects o f a second accelerati on mech- 
anism operating at low energy (JGallant et fflZ.ll2002l ). such as 
cyclotron ac celeration stimulated by a heavy ion compon ent 
of the wind JHoshino et al.lll992l : lAmato fc Aronsll2006[ ). or 
more generally any component with Lorentz factor greater 
than 7uj . We use the spectral continuity to argue that such a 
two component model is unlikely, implying some additional 
piece of physics needs to be added to the shock acceleration 
model to account for these broken power-laws, in addition 
to the extra physics needed to account for the total injection 
rate. 



3 NEBULAR EVOLUTION 

There are three phases of interest in the evolution of a PWN. 
In this section we briefly describe them and explain how we 
model each of them. 

We consider objects of relatively young age (less than 
40,000 yr) with the PSR still inside the SNR, and assume 
sphe rical symmetry, negle c ting t he pulsar proper motion 
(see Ivan der Swaluw et al.l (|200J) for a discussion of this 
point). We also assume, in deriving the PWN evolution, 
that radiation losses of the particles are negligible. Given 
the weak dependence of the PWN radius on the PSR lu- 
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minosity, and the fact that radiation losses only affect the 
particle energy content and not the magnetic energy content, 
we do expect this to be a good approximation. This allows 
us to use analytical formulae, and substantially reduce the 
computational requirements. 

3.1 SNR evolution 

In order to model the evolution of a PWN inside a SNR 
we first need to solve for the evolution of the SNR it- 
self. An excellent semi-analytic model for the evolution 
of an SNR withou t an e mbedded PWN was provided by 
iTruelove fc McKed (jl999l l. Using pressure balance between 
the shocked ISM downstream of the outer forward shock, 
and the shocked ejecta beyond the reverse shock, they solve 
for the evolution of both the outer forward shock and 
the in ner reverse shock from th e early ejecta dominated 
phase (|Hamilton fcSarazinI 119841 ) to the later Sedov phase 
ijOstriker fc McKedl 19881 '). It is shown that the evolution can 
be cast in dimensionless form by using the following charac- 
teristic variables: 



r> »#l/3 -1/3 

Rch = MJ- po 



tch — -C-.q/vAJp,- Po 



^1/2 ,,5/6^ 



(2) 
(3) 



where po is the density of the ISM (assumed to be uni- 
form), EsN is the kinetic energy of the supernova explo- 
sion, whose canonical value is 10^^ erg, and Mej is the mass 
of the ejecta. The only free parameter in the Truelove & 
McKee model is the mass distribution in the ejecta. Unfor- 
tunately, this is not directly constrained by observations, 
and choices in the literature are mostly based on theoreti- 
cal assumptions: for example IChevalierl (|l982r ) suggested a 
distri bution with an inne r plate au and outer steep profile; 
while ITruelove fc McKed (|l999r ) focus on the case of self- 
similar ejecta with a density distribution pej x r~°'t°'~^ 
and a velocity profile v^j oc r. Many recent numerical 
works on young PWNe have indeed adopted the latter 
distribu tion, with great success in r eproducing the obser- 
vations llvan der Swaluw et al. 2001 ;_ van der Swaluwl 20031 : 



iBucciantini et al.ll2003l . |2004| : iDel Zanna et al.ll2004l . |2006| ). 
Once the evolution of the SNR is known, it can be used 
to constrain the evolution of the PWN. 

3.2 Free expansion phase 

The first phase of the PWN evolution is generally referred to 
as free- expansion phase, and has the PWN expanding inside 
the cold ejecta. This phase lasts for a few thousand years, 
until the PWN reaches the reverse shock. In this phase, 
the evolution of the PWN is independent of the evolution 
of the SNR shell, because no contact has been established 
between the two yet. In the case of self-similar ejecta and 
constant PSR luminosity L a solution for the radius of the 
PWN as a function of time has been known for a long ti me 
jChevalier fc Franssoiilll992l : Ivan der Swaluw et aLlbOOll ): 



i?(t) 



1 3 + Q 



M-y't^- 



(4) 



In the more general case in which the PSR luminosity 
changes in time according to 

m^Lo/il+t/rf , (5) 



IBucciantini et al.l l|2004r ) have shown that it is possible to 
find an analytic solution and they provide a series expansion 
for it, showing that already the first few leading terms of the 
series provide an excellent approximation. In this work we 
adopt this approach, that allows us to properly include in 
our model the spin-down properties of the PSR. 

The free expansion phase lasts as long as the radius of 
the PWN is smaller than the radius of the reverse shock com- 
puted using the Truelove fc McKee model. Once the PWN 
reaches the reverse shock, the reverberation phase begins. 



3.3 Reverberation phase 

Once the PWN has reached the SNR reverse shock, it 
is in contact with the shock heated ejecta. The pressure 
of the hot ejecta is higher than the pressure of the rela- 
tivistic material inside the PWN. As a consequence, the 
PWN expansion halts, and the system starts contracting 
(|van der Swaluw et"ai1l200ll '). 

An analytic model for this phase is not available, how- 
ever the evolution ca n be treated wi thin the so called thin- 
shell approximation (|Giulianil Il982l ): the evolution of the 
PWN outer radius is described in terms of the evolution 
of a thin shell of material enclosing a mass Msw equal to the 
mass of ejecta that has been swept up by the PWN up to the 
beginning of the reverberation phase. This shell is bounded 
on the inner side by a hot relativistic plasma with pressure 
Pin, and on the outer side by the shock heated ejecta with 
pressure Pout: 



Ms^R{t) = A-KR[tf[P^„[t) - Pout{t)) . 



(6) 



The value of the PWN pressure is computed numerically 
using energy conservation: 



-finvCJ — ^inyj^r) T 



47ri?(i)4 



L(t)R{t)dt 



(7) 



where tr is the time at which the reverberation phase begins. 
It is less clear what Pout should be. A lower limit is given 
by the pressure that the shock heated ejecta would have in 
the absence of a PWN. An upper limit corresponds to the 
pressure of the Sedov solution. Numerical simulations show 
that the interaction of the SNR with the PWN leads to 
additional heating of the ejecta, due to sound waves that are 
launched inside the SNR du r ing th e reverberation phase (see 
Fig. 3 of IBucciantini et al] l|2003l )'). On the other hand, it 
takes longer for the SNR to relax to the Sedov solution, while 
during the first compression the value of Pout is found to be 
close to about 50% of the Sedov value (|Bucciantini et al.l 
l2003h . In our model we use this fiducial fudge facto r. 

In ID simulations (|van der Swaluw et al.ll2001l ) the re- 
verberation phase is characterized by a series of compres- 
sions and expansions of the nebula, until the system relaxes 
to the Sedov- Taylor phase, which is finally established once 
the PWN reaches pressure equilibrium with the SNR. This 
behavior, which resembles a damped oscillator, is an arti- 
fact of t he ID geometry. In m ore realistic multidimensional 
regime l|Blondin et al.1 l200ll ) , the evolution of the PWN 
during the reverberation phase is subject to strong Rayleigh- 
Taylor (RT) instabilities and efficient mixing of the relativis- 
tic material with the SNR matter. From a dynamical point 
of view this mixing acts as a viscous term on the evolution of 
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the nebula, and one might expect, instead of a series of oscil- 
lations, an almost complete relaxation to the Sedov- Taylor 
solution after the first compression. 

Of course in the presence of mixing and the related 
dumpiness, the volume occupied by the relativistic plasma 
will not be directly related to the radial extent of the neb- 
ula. However Eq. Qwill still hold if one interprets R{t) as an 
effective radius related to the total volume of the relativistic 
plasma, 4TvR{t)^ /3, rather than as the radial extent of the 
nebula (which in general might be larger because of dumpi- 
ness). In 1-zone models this effective volume is also all that 
matters for computing particle adiabatic losses. 

In our model the reverberation phase ends after the 
first compression, once the PWN pressure reaches the value 
of the Sedov-solution for the SNR. 



3.4 Sedov- Taylor phase 

Once the pressure inside the PWN reaches the value proper 
of the Sedov solution corresponding to the SNR forward 
shock, the Sedov- Taylor phase begins. This usually hap- 
pens at an age ~ 10'' yr. By this time, because of spin- 
down, the PSR luminosity can in general be neglected in 
the evolution of the system. The speed of the forward shock, 
Vfa, and the post- s hock pressure are both given by the 
iTruelove fc McKed l|l999r ) model. The pressure inside the 
SNR is then assumed to be equal to the central pressure of 
the Sedov solution Pout ~ 0.5po«/s- 
The nebula evolves according to 

R{tsfPn{ts) = R{t)*P,n{t) = R{t)*Pout{t) , (8) 

where ta is the time at which the Sedov- Taylor phase begins 
and, again, the radius R{t) is related to the total volume of 
the relativistic plasma rather than to the radial extent of 
the nebula. 

One last comment is in order about the role of RT in- 
stabilities. While during the reverberation and Sedov-Taylor 
phase, re-interpretation of the radius appearing in Eqs. |6l 
[7] and [8] in terms of effective volume is necessary, during 
the initial free-expansion phase, RT and mixing are not as 
important: iBucciantini et al.l l|2004bl ) and Ijunl 1 19981 ) have 
shown that, in general, the radius derived from the ID solu- 
tion, is at most 10-15% smaller than the true radial extent 
of the nebula. We might then conclude that at least for very 
young objects the ID radial model provides a reliable esti- 
mate both for the radial extent of the nebula and for the 
volume occupied by the relativistic fluid. 



4 PARTICLE EVOLUTION 

Once the evolution of the nebula is known one can compute 
the evolution of the particle distribution function. 
The energy of a particle is evolved according to 



dEjt) R{t) 



dt 



E(t) 



4crf 



E\t) 



B{tf 
87r 



+ U(t) 



(9) 



where B{t) and U{t) are the magnetic field and the back- 
ground photon energy density in the nebula, respectively, 
ar is the Thompson cross section, m the particle mass. The 
magnetic field in the nebula is computed assuming that the 
ratio of magnetic to total energy in the nebula tjm is constant 



in time (0 < r)M < 1), while the total energy is computed 
considering the pulsar injection and the adiabatic losses of 
the ID model as described in Sec. |3] The background pho- 
ton energy density includes different contributions: CMB 
and starlight, which are constant in time, and synchrotron 
(SYN) emission, which is computed, instead, together with 
the nebular evolution. We have verified that Inverse Comp- 
ton (IC) losses are in general negligible with respect to SYN 
losses, and SYN-IC is important only in young compact ob- 
jects, where the magnetic field is stronger, and the SYN 
emissivity higher. 

Given a particle injected at time to with energy Eo, one 
can solve Eq.|9]to derive the energy E{t, to, Eo) of the parti- 
cle at time t. An analytic solution of Eq. [9] is not in general 
available, neither is its inverse, and the energy evolution has 
to be computed numerically. Once the evolution of the en- 
ergy is known, then it is possible to compute the evolution 
of the particle spectrum according to 



NiE,t) 



' ■ Bt 

N{Eo,to)-^{E,Eo,t)dEo , 
oE 



(10) 



where N{E) is the particle injection rate per unit energy 
interval. 

Since no analytic solution is available for Eq. [O] also 
this equation must be solved numerically. Given the very 
short synchrotron lifetime of the particles at the high energy 
end of the spectrum we use a Lagrangian scheme in energy 
space. The energy space is originally divided into energy bins 
whose extremes are given by Ei{t). These evolve according 
to Eq. [9] When a bin moves to energies lower than 50 keV, 
it is removed from the sample. New bins are continuously 
added at the high energy end of the distribution function, as 
the old ones evolve to lower energies. The time step is then 
dictated by the requirement of sufficient energy resolution 
in the high energy portion of the spectrum. 

The total number of particles inside the i-th bin, at each 
time, is then computed by integrating numerically the fol- 
lowing equation: 

N,{E,+r;2{t),t)^ 
Ni{E,+y2{to),to)+ dt' / dE N{E,t') , (11) 

Jto J Ei(t') 

where to is the initial time at which injection in the bin 
begins. This approach guarantees conservation of particle 
number. 

In principle the functional form of N{E, t) can be arbi- 
trarily chosen. 

For the the pairs, we assume a broken power-law. Given 
that injection is due to the PSR wind, and that the only 
energy scale available in the wind is the PSR voltage, e„ = 
e$, we have decided to use a scale-free model for a broken 
power-law: 

iV(£,f)-Co(i)| (^/^^)-.. fore,<£<6„ (^2) 

where 71 < 2 < 72, tc is the peak of the injected energy dis- 
tribution EN{E, t) and tm is the low-energy cutoff, usually 
found to be small compared to Ec. The scale free approxi- 
mation implies that the ratios u^ — em/e„ and fie = ^c/^v, 
axe constant in time, /ie and Ve are free parameters of the 
model, and their value is obtained from a fit to the data. 
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Once /Xe and v^ are fixed, Co{t) is given by energy con- 
servation. We have: 



r;,L(f) = / N{E,t)EdEi 



72 —71 



(72-2)(2-7i) 



Co{t)el (13) 



with r^e the fraction of the PSR luminosity injected into 
pairs. The particle number fiux in the wind is 



Nit) 



N{E,t)dE: 



Co{t)e 



71 



1 Ve™J 



(14) 



and the average energy/particle in the spectrum (|12|) is, for 

rieL{t) 



(E) 



N{t) 
(72 -71) (71 



= JwirieC 



71) V Ec 



(72-2)(2-7i)^^Vecy • ^^^^ 

In expression 1151 we have identified the average energy/pair 
in the spectrum with the upstream flow energy/pair, but we 
have not speciflcally used the shock jump conditions. This is 
because, as we will discuss in iJTl the distribution (I12|l implies 
some additional mechanism for energy redistribution, which 
may or may not be directly connected to the termination 
shock. From Eq.[TS]it is clear how it is possible to estimate 
the PSR multiplicity through a fit of the model to the PWN 
emission. 

If the wind is characterized by a single value of the 
Lorentz factor, then 'y^meC? ~ ec, rieL{t) — me(?'yn]N{t). 
However, the approximation of a single Lorentz factor might 
not be correct, in this case the ratio tm/tc cannot be a- 
priori determined, but is found from requiring a fit to the 
entire spectrum. In particular, once an upper limit for em is 
derived from fitting the low frequency radio emission, this 
can be used to infer the particle number flux and estimate 
the average Lorentz factor. 

We consider the possibility that a minority by number of 
high energy particles in the equatorial return current (either 
positrons, ions or electrons, depending on return current 
sign) are part of the PSR wind and carry a fraction r]p of the 
energy flux. Focusing on the case when ions are the implied 
particles, one might look for signature of their presence. For 
simplicity we consider them injected with a mono-energetic 
spectrum centered at an energy ^^mpC? (equivalent to as- 
suming that they are moving at the same Lorentz factor as 
the pairs injected at ec). Ions are affected by different loss 
mechanisms compared to pairs: while radiation losses are 
negligible, p-p scattering and diffusion outsi de the nebula 
might be important. Following the approach of lAmato et al.l 
1120031 ) we assume that both diffusion and p—p scattering re- 
move particles from the distribution function. Due to these 
effects the number of ions of energy E between time t and 
t + dt changes by: 

N{E{t + dt),t + dt) = iV(£'(t), t)e~*<^''"''^+'/"<*'-f ^ (16) 

where E{t) is given by Eq. [O] and where the characteristic 
p ~ p timescale is 



5{aontc) 



(17) 



with CTo = 5 X 10 ^^ cm^, and nt the target proton density. 
While the typical diffusion timescale is 



TMf 



where tl is the ions' Larmor radius. 



5 EMISSION 

There are three channels of non-thermal emission in PWNe. 
From Radio frequencies to MeV photon energies, the dom- 
inant emission process is synchrotron radiation by the ac- 
celerated electrons and positrons gyrating in the nebular 
magnetic field. At higher energies, there are two possible 
contributions: Inverse Compton (IC) emission by the high 
energy pairs, and gamma-ray emission from neutral pion 
decay, following "p-p" scattering between accelerated and 
target protons. 

The synchrotron luminosity is computed using the 
"monochromatic" approximation. The power per unit fre- 
quency emitted by a single particle with energy 'yrUef? is: 



5(^,7) 



6 on 



(19) 



where ctt is the Thompson cross section, B is the nebular 
magnetic field, and 



0.29^^57" 



(20) 



A-Kmc 
is the characteristic emission frequency. 

Eq. [19] has then to be integrated over the pair distri- 
bution function (see previous section) in order to derive the 
total nebular synchrotron luminosity at a given frequency. 

I C is computed using the full Klein-Nis hina cross sec- 
tion (jJoneslll968l : iBlumenthal fc Gouldlll97d ') and assuming 
scattering over three main target photon flelds: 1) CMB, 
modeled as a prefect black-body, with typical temperature 
2.75 K; 2) synchrotron emission from the PWN, consistently 
computed as described above; 3) "galactic background", 
which in principle includes both galactic starlight and any 
local contribution in the optical-IR band. This third term is 
modeled as a suppressed black-body: both the temperature 
(in the range of a few hundreds/thousands Kelvin) , and the 
suppression factor (of order 10"^"), are in principle free pa- 
rameters that can be adjusted. The role of this contribution 
will be discussed in more detail when dealing with specific 
objects. 

Gamma ra y emission from neu tral pion decay is com- 
puted following lAmato et al.l (|2003h . In the scaling regime a 
mono-energetic proton distribution with energy Ep leads to 
the following pion injection spectrum 



dN^ 
dtdE^ 
where 



— Kt^E^ g-^{ETr/Ep) 



g^{x) = {1 - x)''^ + e-'^yi.34 



(21) 



(22) 



and the constant Kt^ is found through the normalization 
condition: 



mE.){^ 



dN^ 
dtdE, 



■ E-jY cLEty 



(23) 



where N{Ep) is the number of protons with energy Ep. 

The number of photons emitted per unit time and en- 
ergy interval is then 



tlc 






(24) 
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The total pion emissivity is then computed integrating 
over the total proton energy distribution. The density of 
target protons nt that enters Eq. [17] is the main unknown; 
if the protons are well confined within the PWN, then a 
reasonable assumption would be to use the swept-up ejecta 
mass; if they can leak out then one should consider the total 
ejecta mass. A realistic model for proton diffusion in PWN- 
SNR systems is far beyond the scope of this paper: we keep 
nt as a free parameter and discuss the implications of its 
assumed value whenever relevant. 



6 RESULTS 

In this section we show the results obtained when apply- 
ing our model to a few objects of different ages, both young 
and old. For each object, we will not attempt at obtaining 
the best possible fit of the data, which in our view would 
not be very significant given the intrinsic limitations of a 
one-zone model. In fact, our primary interest is to put some 
constraints on the unknown parameters and to clarify which 
observations would be more relevant to improve our under- 
standing of these systems. At the same time we will discuss 
the limitations of this approach and make clear which part 
of our results should be considered as really reliable. 



6.1 The Crab Nebula 

The Crab Nebula is the system for which we have the best 
constraints, both regarding the pulsar injection properties, 
and the luminosity at different wavelengths. A vailab l e data 
extend from low radio frequencies (Baldwin Il97ll: iBaara 





1972), through mm/IR (iMezgCT et al. 1986; Bandiera et al.| 
20021) ■ optical/UV ra diation dVeron-Ccttv fc Woltj er 199^ 
Hennessv et~al] 1 19921 ), X-ravs (JKuiper et all l200ll) and 7- 
rays from MeV to TeV energies ([Aharonian et al.l l2004l : 



[Albert et aP I2OO8I : lAbdo et aL l bOlol ). The source distance 
is es timated to be 2 kpc, and the nebular volume is 30 
pc^ l|Hesterl I2OO8 ) , corresponding to a spherical radius of 
about 6 ly. This is the only system for which we know 
the exact age (950 years). We also have a relia ble estimate 
of the braking index (2.51 implying /3 = 2.33, iLvne et al.l 
1(1993)), which, when combined with the present luminosity 
L ~ 5 X 10^* erg s~^, allows us to derive an initial luminosity 
L(0) = 3.5 X 10'^® erg s~^, and a characteristic spin-down 
timescale r ~ 730 years. 

Much more uncertain are the properties of the SNR 
which should be surrounding the nebula. Several attempts 
at detecting both the forward and the reverse shock have 
only provided upper limits. More information can be de- 
rived by studying the filamentary network surrounding the 
PWN and comparing the results with theoretical expecta- 
tions based on hydrodynamical models of the i nteraction be- 
tween the nebu l a and the surrounding ej ecta IjHester et al.l 
ll996l : [jun|[l99i : lBucciantini et al.ll2004bl ). The visible rem- 
nant contains at least 1-2 Mq of He-rich line-emitting gas 
and the PWN expands at a speed of about 1.5 — 2 x 10^ 



km 



Traditional models for the evolution of the Crab 



Nebu l a llChevalier fc FranssonI 19921: van der Swaluw et al.l 
I2OOII : iBucciantini et al.l I2OO J : iDel Zanna et al.l bOO^ I have 
all assumed a flat density profile of the ejecta. Steeper 



ejecta (a ^ 1 in Table [T| give a mass in the filamen- 
tary network higher than what is observed. Assuming flat 
ejecta and a canonical SN energy of 10^^ erg, in order to 
reproduce the proper size of the nebula a mass Mej = 
9.5M(:) is needed, very close to the most recent estimate by 
iMacAlpine fc Satterfleldl (120081 ). The lack of any clearly de- 
tected SNR shell prevents any reliable estimate of the ISM 
density and we simply assume a flducial value of 0.1 cm~^. 

The PSR and SNR parameters are listed in Table [T] to- 
gether with the parameters that describe the particle energy 
distribution at injection. The low energy spectral index 71 
is fixed according to the radio data, while the high energy 
index 72 is chosen in order to minimize the X-ray residuals. 
We want to stress here that the X-ray emission is concen- 
trated in the central region of the PWN and it is strongly 
affected by the details of t he fiow dy namics just downstream 
of the termination shock (jDel Zann a ct al. 2004, 2006). It is 
unrealistic to expect that a simplified one-zone model can 
provide an accurate description of the high energy spectrum. 
A one-zone model can at most provide an indication of the 
best power-law index that can fit the data, but more real- 
istic multidimensional models are necessary to address the 
emission properties in the X-ray band. Two features are in- 
teresting to notice: 1) 72 is higher than the value 2.23 typical 
of relativistic Fermi shock acceleration with isotropic scat- 
tering in the fiuid frame; 2) there are indications that a single 
power-law cannot reproduce the complete set of X-ray data 
points. 

Several features are present at high energies, but it is 
not clear if they are intrinsic (a sugges tion in this sense has 
been put forward bv lVolpi et al.l (I2OO8), where they present 
a multidimensional model of the emission that shows similar 
features) , if they are simply due to calibration issues among 
different instruments or even for one instruments at different 
energies, or if they are due to tempo ral variability of the 
emission properties (|Bucciantinill2008l ). 

The present model can reproduce the synchrotron part 
of the spectrum within few % accuracy (greater discrep- 
ancies being present at high energies as discussed above). 
At present, our particle distribution function does not in- 
clude any smooth high-energy cut-off, the reason being the 
very short lifetime of high energy particles that prevents ef- 
ficient numerical integration of their evolution. However a 
good match with the data can be obtained by introducing a 
posteriori an exponential cut-off in the particle energy dis- 
tribution. This refiects in a cut-off in the emission spectrum 
at a frequency around ~ 10^^ Hz. Our model requires that 
80% of the pulsar spin-down luminosity goes into acceler- 
ated pairs. This does not leave mu ch energy to be ca rried 
by higher energy particles (i.e. ions IjAmato et al.ll2003r )). In 
principle the presence of ions might have observable conse- 
quences in the TeV band, however even assuming the re- 
maining 20% of the spin-down energy is all carried by ions, 
their contribution to the gamma-rays emission turns out to 
is negligible. On the other hand, if the return current is made 
by a high energy lepton beam (Arons, in preparation), the 
radiative consequences of this latter scenario have not been 
investigated, however, again, they are not expected to give 
any appreciable contribution to gamma-rays 

TeV emission is due to IC scattering by the pairs, for 
which the main target is the nebular synchrotron emission 
in this case. We estimate a magnetic field B ~ 200/iG, 



N. Bucciantini, at al. 



Table 1 . Values of the parameters used in the modehng. 



Parameter 


Symbol 


Crab 


3C58 


B1509-58 


Kes75 


Supernova explosion energy (10^^ erg) 


EsN 


1 


1 


7 


2.1 


Mass of the ejecta {Mq) 


Mej 


9.5 


3.2 


4.0 


16.4 


ISM density (cm-^) 


Po 


0.1 


0.01 


0.001 


2 


Ejecta density index 


a 





1 


1 





Initial pulsar luminosity (lO'^^erg s~^) 


Lo 


35 


0.73 


49 


1.66 


Spin-down time (yr) 


T 


730 


3280 


114 


226 


Braking index 


P 


2.33 


2 


2.087 


2.12 


Age (yr) 


t 


950 


2100 


1570 


650 


Fraction of L that goes into pairs 


Ve 


0.8 


0.75 


0.7 


0.95 


Low-energy injection index 


71 


1.5 


1.2 


1.2 


1.7 


High-energy injection index 


72 


2.35 


2.82 


2.14 


2.3 


Peak energy ratio 


/ie 


1.54x10-5 


6.25xl0-'3 


3.33xl0-« 


10-4 


Minimum energy ratio 


Ve 


1.1x10"^ 


8.7x10"* 


2.8x10-8 


2 X 10-^ 


Break energy (eV) 


f-c 


4x10" 


5x10" 


1.5x10" 


4x10" 


Fraction of magnetic energy 


r\M 


0.11 


0.5 


0.53 


0.005 


Fraction of L that goes into ions 


Vp 















Luminosity 




10'= 10"" 10" 

Frequency (Hz) 

Figure 1. The C rab Nebu l a integ r ated em i ssion sp ectrum. Data 
points are from iBaldwinI lll97]); iBaarj lll972l'): iMezger et aO 
ill986l'): ISandieraet al.l (|2002|'I: Iveron-Cettv fc Woltied lll993l): 
iHennessv et al.l lll99lV: iKuiper et al.l Il200lh: lAharonian et all 



Ii2004 ): lAlbert et al.l bOOSi) : lAbdo et ahl bOlOl ') Ctriangles are 
EGRET points, squares are Fermi points). Solid line is the total 
luminosity. Dashed lino is the IC-CMB, dotted line if the IC-SYN. 



slightly l ower than previous estimates (|Aharonian et al.l 
I2OOJ; iHestcr 2008). Overall the broad-band spectrum is very 
well reproduced, with the largest discrepancies limited to the 
EGRET data points around lOMeV. These points have large 
uncertainties and are likely affected by calibra tion issues ( 
("Abdo et al."2009')). Indeed recent Fermi data (lAbdo et al.l 
L2OIO) arc consistent with our model curve. 

iBaldwinI l)l971f ) has shown that the radio power-law 
spectrum of the Crab Nebula extends down to the iono- 
spheric cut-off at ~ 30 Mhz. One can use this piece of in- 
formation to derive an upper limit on the minimum particle 
energy at injection. In order for the radio spectrum to ex- 
tend to those energies as an uninterrupted power-law, we 
need to assume £m/ec < 7 x IQ-*. From Eq. [T5]we then find 
for the wind Lorentz factor 7^, < 5 x 10*. This has to be 
compared with the typical Lorentz factor of the particles at 



€c, which is 7.4 x 10^, and with the minimum Lorentz factor 
of injected particles that is < 500. The estimated value of 
7u, translates into a lower limit on the pair multiplicity of 
K > lo**. 

Our model allows us to compute a posteriori the energy 
radiation losses. We do this in order to verify to what degree 
the adiabatic approximation for the evolution of the PWN 
radius is correct. For the Crab Nebula we find that about 
half of the total energy injected into the nebula in the form 
of pairs has been lost via synchrotron emission. However the 
dependence of the radius on the pulsar luminosity L{t) is in 
general weak (scales as L{tY'^ for constant luminosity), so 
we expect at most modification of order 20% in the radial 
evolution, and maybe in the adiabatic losses. Such value is 
well within the simplifications and the approximations of 
the model. In order to properly take into account radiation 
losses, one needs to abandon any analytic solution for the 
dynamics, and solve the coupled system of equations for the 
dynamics and the emission simultaneously. 



6.2 3C58 

3C58 is another example of a young PWN, which shares 
many similarities with the Crab Nebula. It has a typi- 
cal non-thermal spectrum extending from Radio to X-rays 
llSalter et al. 1989: Torn et al.ll2000l : K?reen fc Scheueilll992l : 
ISlane et al.. .2004 . .2008. ) . It shows clear evidence of energy 
injection from the cent ral pulsar PSR J0 205+64 in the form 
oi a jet-torus structure l|Slane et al.ll2002l ). Recent SPITZER 
measurements of the PWN IR luminosity have showed clear 
eviden ce for a possible injection break, like in the case of 
Crab (jSlane et al.ll2008l ). At a typical distance of 3.2 kpc 
(IChevalierl l2005h , its size is about 5x9 pc, equivalent to 



a volume of ~ 140 pc (corresponding to a spherical ra- 
dius ~ 3. 3 pc) . An association with the SN 1181 has been 



proposed l|Stephenson fc Greerj|2002l ). However recent mea- 
surements of the exp ansion of the filamentary structure sur - 
rounding the source l|Bietenhol j|2006l : iRudie fc Fesenll2007l ) . 
interpreted through a simple model of the PWN-SNR sys- 
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Table 2. Values of the parameters corresponding to different 
models for 3C58 



3C58 Luminosity 



/3 


a 


t(yi-) 


M 


y (Mq) 


VSH 


(kms-1) 


2 





2000 




5.4 




270 


2 


1 


2100 




3.2 




305 


2 


2 


2250 




1.2 




355 


2.33 





2050 




5.5 




270 


2.33 


1 


2180 




3.3 




307 


2.33 


2 


2350 




1.3 




360 



tern l|van der Swaluwll2003l : IChevalieiil2005l ). have ruled out 
this possibility. 

The pulsar present spin-down power is L ~ 2.7 x 10^^ 
ers s~^, with a ch aracteristic dipole age fc ~ 5390 yr 



JMurrav et al.ll2002l ). Neither the true age nor the PSR brak- 
ing index are known, and this prevents any reliable determi- 
nation of the initial spin-down power. However it is possible 
to obtain some constraint from the information available on 
the exp ansion of the PWN a nd surrounding SNR. In par- 
ticular, iBocchino et al.l (|200l|) measured a mass in the fila- 
ments of ~ O.liMQ, and thermal emission corresponding to 
a typical shock velocity ~ 340 km s~^. Assuming that the 
filamentary structure corresponds to the swept-up shell of 
ejecta (in analogy with the case of the Crab Nebula) then 
we can use our model for the system evolution and obtain 
an estimate of both the system age and the PSR braking 
index depending on the density profile of the ejecta. 

In Table (2] we present a set of models derived for dif- 
ferent values of the unknown braking index (3 and 2.5), and 
for different density profiles of the ejecta. In all the dif- 
ferent models the swept-up mass and the nebular volume 
are consistent with extant measures. The shock speed varies 
slightly, but is in all cases consistent with measurements, 
within the approximation of a spherical model. We again 
assume a canonical SN explosion energy of 10^^ erg. As in 
the case of Crab the lack of any clearly detected SNR shell 
prevents any reliable estimate of the ISiM density and we 
simply assumed a fiducial value of 0.1 cm~^. 

By looking at Table [2] one immediately realizes that 
the inferred age is not very sensitive to either the braking 
index or the assumed ejecta structure. The reason for this is 
that in all cases the age is much less than the characteristic 
age, and the pulsar spin-down power in the past was not 
much larger than the present value. It is interesting to notice 
that in all models the shock speed tends to be higher for 
steeper ejecta, favoring the case r~^ . However this would 
result in a very small mass of the ejecta, inconsistent both 
with estimates based on the observed filamentary knots (~ 8 
Mq) (jRu die & Fescn 2007), and with stellar evolutionary 
models. 

It is our opinion that cases with a density profile which 
is either flat or r~^ are more hkely, even if they are associ- 
ated with a lower shock speed: this is because steeper ejecta 
lead to unreasonably small ejecta masses. In all cases the 
age estimated for the system rules out any association with 
SN 1181. However the filamentary structure, that we have 
interpreted as mostly due to swept-up ejecta . is rather com- 
plex and with various velocity components ( Rudie fc FesenI 



10^ 
10^ 



10= 




Frequency (Hz) 



Figur e 2. 3C58. Data points are: radio fro m from ISalter et al.l 
(1989'! iGreen fc ScheueJ lll992h . X-ray from ISlane et al.l l|20ojr 
IR from ISlane et al.l ll2008h (diamonds are IRAS and SPITZER 
points). Solid line is the total luminosity. Dashed line is the IC- 
CMB, dotted line is the IC-SYN. 



[20071), so that it cannot be ruled out that part of this fila- 
mentary network might be associated with the 1181 event. 

The question that naturally arises is whether spectral 
fitting allows to discriminate between the models listed in 
Table [21 Unfortunately the answer is negative: fits to the 
observations lead to similar results in terms of the required 
injection parameters for all of these models. The reason is 
that all of them correspond to similar underlying ages, sim- 
ilar integrated spin-down power (i.e. total injected energy) 
and the same size (same adiabatic losses). The parameters 
reported in Table [T] correspond to the case of ejecta with 
a = 1 and braking index 3. Spectral fitting is shown in 
Fig. [2] Very similar values of the parameters are found also 
for a braking index of 2.5 (/3 = 2.33). 

If we compare the particle spectrum at injection with 
that of Crab, we find that the spectral index is now flatter at 
low energies and steeper at high energies. Again the required 
pair injection efficiency is very high ~ 80% and the ratio 
fie ~ 8.3 X 10~^ is only a factor 2 lower than for Crab. The 
inferred magnetic field is 4 3AtG, higher than other recent es- 
timates (jSlane et al.ll2008h . The ratio between the magnetic 
energy and the total energy is ~ 0.5, this is quite high and 
leads to a value for i]e + tjm > 1 which apparently violates 
the total energetics. The discrepancy is however ~ 30% and 
we consider it to be within the approximations of our model 
(particularly critical in this respect are the simplifications 
related to the assumption of spherical geometry) and the 
uncertainties of the data (especially the estimated distance, 
and hence volume). Models with rje + tjm = 1 always un- 
derpredict the radio emission, independently of the value 
chosen for rjM- 

Within the present model, we consider our estimates of 
rje and r]M rather solid: attempts at fitting the data assum- 
ing a lower magnetization and a higher particles' content 
result in underproduction of the radio emission, while larger 
magnetization leads to the break frequency moving below 
10^^ Hz, in contrast with observations. One final remark we 
want to make about tjm is that the idea of a relatively high 
magnetization in this PWN is also supported by its sub- 
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stantial elongat ion, created by hoop stresses of the toroidal 
magnetic field (|Begelman fc Lilll992l ). 

The existence of radio measurements down to 100 MHz 
allows us to derive also in this case a lower limit on the pair 
multiplicity: we find that 711, < 3 x 10* to be compared with 
the typical Lorentz factor of the particles at tc, which is 

1.2 X 10^, and with the minimum Lorentz factor of injected 
particles that is < 1600. This implies k > 5 x lO''. AU the 
latter values are very close to the corresponding estimates 
for Crab. 

6.3 MSH 15-52 

The system formed by the pulsar PSR B 1509-58 and the 
SNR MSH 15-52 is striking for its relatively large size com- 
pared with its pres umed age. Its d istance is estimated to be 
5.2 ± 1.4 kpc (jGaensler et al.lll99& '). Assuming the central 
value as a fiducial distance (which we will do in the follow- 
ing), the SNR radius turns out to be ~ 83 x 10^* cm. PSR 
B 1509-58 is a young pulsar. Its current spin-down power is 
1.8 X 10^^ erg s~'^, with a characteristic dip ole age of 1550 
yr (JKaspi et a.l.lll994l : lLivingstone et al.ll2006l ). This is one of 
the very few pulsars for which the braking index is known: 
P — 2.087. Unfortunately, our ignorance of the true age of 
the system still prevents knowledge of the initial spin-down 
power, Lo'. this can only be determined as a function of the 
assumed age, t. A maximum possible age of the system can 
be estimated from the braking index and the characteristic 
dipole age and the result is t < 1690 yr. For such an age, 
the size o f the SNR implies a high Esn/M^j ratio, typical of 



Table 3. Different age rages for models for B1509-58 



SN Ib/c (tCaensler et al.lll999l:lTsvetko v"2002': ' Mazzah et al 
2OOOI: llwamoto et al.l2000l : lMazzah et al...200a : ITanaka et al 



20091 ). The size of the PWN is poorly constrained because 
of the high radio foreground from the SNR. Spectr al model- 
ing ( Nakamori ot al. 2008; Chevalier 2005; Du Ple ssis et al.1 
[1993), however, suggests that the cooling frequency lies just 
below the CHANDRA band (see also the present discus- 
sion), which is also in agreement with the size of the neb- 
ula at TeV energies being comparable with the CHANDRA 
size. If cooling starts being important at frequencies no lower 
than the X-rays, then one can estimate the size of the PWN 
from CHANDRA images: the result is a PWN radius of 
~ 17 X 10^* cm, for the assumed distance. 

Knowing the braking index, and assuming an age (i) 
and an ejecta density profile (a), we can use our model to 
determine which values of Esn and M^j return the assumed 
radius for both the PWN and SNR. The only other param- 
eter that enters this calculation is the local density of the 
ISM. We assume po = 0.001 cm~3: a higher density would 
hint at even higher values of Esn/Msj and in any case does 
not allow to match both radii; a lower value marginally af- 
fects the results because the system turns out to be in the 
ejecta dominated phase anyway. 

If one accepts the hypothesis that B1509-58 was born 
in a SN Ib/c explosion, then one can further constrain the 
model requiring that the mass of the ejecta is in the range 

4 — IOM0, and the kinetic energy released is in the range 

5 — 20x 10^^ erg. In Tab.[3]we report the different age ranges. 
Flat ejecta tend to favor slightly younger systems. In spite of 
the relatively small differences between the ages reported in 
Tab.|3l spectral modeling can prove a useful tool to discrimi- 
nate between different scenarios for this source. In the case of 
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1350-1500 
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1.2 x 10-"^ 
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3C58 the inferred ages for all models were much smaller than 
the characteristic dipole age, so that the injection properties 
of the PSR were almost the same, and it was not possible to 
discriminate using the SED. On the contrary, for B1509, the 
inferred ages are close to the characteristic dipole age, and 
different models correspond to different spin-down histories, 
and different integrated injection energies. In such a situa- 
tion the SED, and in particular the radio data, can be used 
to rule out some of the models allowed by t he dynamics. 

Observations cover the radio GHz band (|Gaensler et al.l 

[1999I) (although the quality of the data in this band is 
not very high) X-rav s from CHANDRA to INTEGRAL 
llGaensler et al.l 12003 : iDeLanev et"ai] I2OO6I: iForot et al.1 
l2006l ). Fermi in the MeV-GeV rang e JAbdo et al."2010bh 
Cangaroo-II and HES S in the TeV (JNakamori e t al. 20oi; 
lAharonian et al.ll2005l ). Only models with Lo > 10^'' erg 
s~^ can reproduce the radio, while lower values of Lo do 
not satisfy the energetic requirements. Tab. [3] shows the 
limits imposed on the model by SED fitting. Flat ejecta 
correspond to a very narrow range of possible ages and 
a value of fie ~ 1.2 x 10~® smaller than what is found 
both for the Crab Nebula and 3C58. Vice versa, ejecta with 
Q = 1 correspond to spin-down properties and multiplic- 
ity that are closer to what is found in Crab and 3C58. 
Compared with these two sources, we find a slightly harder 
high energy spectrum at injection 72 ~ 2.1 — 2.16. How- 
ever, this is consistent with the hard spectrum obs erved 
in the MeV-GeV range by Fermi (|Abdo et al.ll2010bl ). In- 
deed the photon index derived from X- ray observations 
(G aensler et al.ll2002l : [beLanev et al.ll2006l ) of the inner re- 
gion is ~ 1.4 — 1.6, and is much smaller than the previously 
used value ~ 2.7 which seems inconsistent w ith both obser- 
vations and energetics (jNakamori et al.ll2008l ). Moreover the 
injection break is located between 5 and 50 GHZ, at a signif- 
icantly l ower frequency than w hat was assumed in previous 
models JNakamori et akllJOol ). 

Steeper ejecta give larger values of fie, as reported in 
Tab. O All possible models require a high acceleration effi- 
ciency, rfe ~ 0.7 and a strong magnetic field, rfM ~ 0.53 im- 
plying B > 20/xG. Again, as in the case of 3C58, there is an 
energetic problem, but as in the previous case, this system is 
strongly magnetized , as also suggested by the presence of a 
strong X-ray jet (see lDel Zanna ot al. ( 2004) for a discussion 
of the dependence of the jet strength on the magnetization). 

The lowest frequency radio point allows us to constrain 
the value of Ve, as it was done in Crab and 3C58. We find 
that even for a = 1, the average wind Lorentz factor is 
jw < 10'', quite small compared with other systems. This 
translates into a lower limit on multiplicity, that has to be 
greater than ~ 2 — 3 x 10^ . 

Fig. |3] shows the spectrum derived from the model pre- 
sented in Tab. [T] From spectral fitting it is clear that the 
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Lorentz factor corresponding to Ec has to be smaler than 
5 X 10* (/ie has to be smaUer than 4 x 10""), otherwise it 
is not possible to reproduce the correct slope and intensity 
in the radio band (the injection break shifts below IGHz). 
As a result rjM must be larger than 0.3, and the magnetic 
field cannot be smaller than 15/xG. Only by measuring the 
PWN emission in the range 10^" — 10^^ Hz, it would be pos- 
sible to discriminate among various models and put better 
constraints on the multiplicity. 

As it has already been noted (jAharonian et al.ll2005l ) it 
is not possible to reproduce the observed gamma-ray emis- 
sion in the TeV range, either with IC on the CMB or the 
average galactic background. The contribution from IC-SYN 
is negligible. Various models have been presented to account 
for such discrepancy. 

One suggestion is that of a lower value of the magnetic 
field in the nebula, which would lead to infer a higher particle 
conte nt, and consequently en hance the gamma over X-ray 
ratio (jAharonian et al.ll2005l ). However, within this model 
this does not seem to be a viable solution: in order to in- 
crease the particle content up to the value required to fit the 
gamma-rays as IC on the standard galactic background we 
would have to violate the energetic fixed by the spin-down 
history by a factor ~ 5. 

Another proposal is that of a possible contribution to 
the TeV fiux of n° decay, in the presence of relativistic 
hadrons (jNakamori et al.ll2008l ). It was however immediately 
reaUzed that this would require so much energy to be put 
in protons as to require millisecond magnetar conditions 
at birth for the pulsar, which seems unlikely judging from 
the present PSR-SNR properties. In addition, the swept-up 
ejecta do not provide a sufficient target density of thermal 
protons. 

The third possibility is that the local photon back- 
ground, in particular the IR, could be much higher than the 
average galactic val ue. Indeed the SNR itself could be the 
origin of this excess (jNakamori et al.ll2008l : IDu Plessis et al.l 
[l995|). The fit presented in Fig. [3] assumes for the local IR 
background a black body with a temperature ~ 400K sup- 
pressed by a factor 3 x 10~^, which corresponds to an energy 
density about 5-7 times higher than the average galactic 
ba ckground . It is interesting, in this regard, the suggestion 
by iHelfandl (|2007l 'l that the correlation between SNR de- 
tected by HESS and HII regions might be due to enhance- 
ments in the local photon background. 



6.4 Kes 75 

The supernova remnant Kes 75 (G29.7-0.3) is a shell type 
remnant which hosts at its center a young PWN powered 
by pulsar PSR J1846-0258. This is a 324 ms pulsar, with a 
characteristic dipole age of 723 yr, a spin-down luminosity 
of 8.3 x 10'^'' erg s~^, and a measured value of the braking 
index /3 = 2.65, which p rovides an upper limit fo r the age of 
the system of ~ 880 yr ([Livingstone et aLll2006f ). 

The estimated distance of the system suffers from large 
uncertainti es. The first estimates pointed to a distance of 
19-21 kpc (jBecker fc He lfand 1984)- However it was imme- 
diately realized that such a large value impli ed very pecu- 
liar conditions for the Su pernova explosion (|Morton et al.l 
l2007l : iHelfand et al] |2003| '1 and the ISM, together with an 
extremely high synchrotron efficiency of the PWN. Recent 
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Figure 3. B1509-58; M SH 15-52. Data p oints are: ra- 
dio f ro m iGaen sIer et al.l lll999|). X-rav f rom i GaensIer et al.l 
l!2002l) IP eLanev e~ ||2006|') iForot et all (120061), 7-ray from 
Nakamori et al.l ll20oi ): lAharonian et al.l l|2005l '): iKuiper et~all 
(il999i) : lAbdo et al.l ( 12010^ ) fdiamonds are EGRET points, 
squares are Fermi points, triangles are COMPTEL/ HESS 
points). Solid line is the total luminosity. Dashed line is the IC- 
CMB, dash-dotted line is the IC on the enhanced local IR back- 
ground at 400K. 



studies have revised the distance to a much smaller value, 
in th e range 5-7 kpc (jEisenhauer et al.ll2005l : iLeahv fc TianI 
120081 ). We will assume a distance of 6 kpc in all our models. 

Assuming a distance of 6 kpc, the radii of the PWN and 
of the SNR are, respectively, R = 1.8 and Rsnr = 9.24 ly. 
The velocity of the fo rward SN shock is Vsn = 3700 km s"'^ 
(jHelfand et al.l I2OO6I ) as estimated via the Si X-ray line. 
This places an upper limit on the age of the system t < 
Rsnr/Vsn ~ 800 yr. Various attempts have been made to 
measure the ejecta mass, and the SN energy; however the 
bulk of the SNR emission comes mainly from two regions of 
the SNR shell, and this tends to bias all estimates. In gen- 
eral, very high masses are inferred both for the s wept-up ISM 
and for the SN ejecta (JBecker fc Helfand 1984 : Morton et al.1 
I2OO7I : iLeahv fc TianI I2OO8I ). hard to reconcile with theoret- 
ical expectations for a SN Ib/c and standard evolutionary 
models. It must be stressed, however, that mass estimates 
are very sensitive to the assumed volume of the emission 
region, which is not well known. 

Proceeding as for previous cases, we will first attempt to 
constrain the possible SNR-PWN parameters through mod- 
els for the dynamical evolution of the system. Given an as- 
sumed age t and an ejecta profile a there is only one set 
of values of Esn, Mej and po that gives the correct size 
of both the PWN and the SNR and the correct forward 
shock speed. Further constraints on the parameter space 
come from the properties of SN Ib/c which limit Mej in 
the range 5Mq < M^j < 16Mq. Cases with a > 1 are ac- 
ceptable only for a system older than 700 yr, and are in 
general unlikely because they require a very low SN energy, 
_Esjv < 0.3 X 10^^ erg. Cases with flat ejecta are admissible 
for a larger range in age from 450 yr to 650 yr. The typi- 
cal SN energy does not seem to depend on the age and is 
Esn ~ 2 X 10^^ erg; younger systems require lower values 
of Mej, and larger values of po, as shown in Tab. |4| Cases 
with marginally steep ejecta, a = 0.5 suggest a larger age. 
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Table 4. Different age ranges for models for Kes75 
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a less energetic explosion, Esn ~ 0.6 x 10 erg, and lower 
ISM density. 

The question arises again of whether it is possible to 
discriminate among the different parameter sets for the dy- 
namics based on the SED. Unfort unately this is n ot the case: 
as it was already realized earlier l|Chevalierll2004 ) . Kes 75 is 
a particle dominated system by large (more than 90% of 
the pulsar spin-down power seems to be converted into ac- 
celerated particles, see below), so the differences in particle 
content due to the different energetics of the various models 
can be easily compensated by small changes in the magne- 
tization to give the same synchrotron emission. 

However, fitting the model to the multi-band emis- 
sion spectrum is non-trivial. The radio emission spec- 
trum implies a low energy particle spec t ral index at injec- 
tion 71 = 1.7 llBecker fc Helfandl Il984l : ISalter et al.l 1 19891 : 
iBock &: Gaenslerll2005l ). The average photon index in the 
CHANDRA band is 1.9 while deep X-ray images have shown 
that the photon index in the vicinity of the pulsar is '--^ 1.8 
dBlanton fc Heffandl[l99^ : ICollins et aLll2002l : iHelfand et all 



l2003l : lNget al.ll200a ) 



Evidence for a p ossible spectral break below 100 
GHz was presented bv iBock fc Gaensleri l2005|), but these 
data seem to be i nconsistent with previous measurements 
l|Salter et al.lll989l ') and a single power-law cannot be ruled 
out. 

Indeed in our model an injection break below lOOGHz 
would require a very hard high energy injection, inconsis- 
tent with spectral information at X-ray frequencies. IR data 
do not provide good constraints IjMorton et al.l 120071 ). IN- 
TEGRAL data above 15 keV show a particularly hard spec- 
trum, not fully consiste nt with an expecte d cooling break 
in the CHANDRA band l|Terrier et al.ll2004r ). but the pulsar 
might contribute to the fiux above 40 keV and be responsible 
for the excess emission at high energies. 

If the X-ray emission in the CHANDRA band corre- 
sponds to freshly injected high energy particles, one must 
infer 72 > 2.6. However we found that it is not possible to 
fit the overall SED when adopting such a soft high energy 
injection. The reason is simple: given that the average spec- 
tral index is close to the one measured at the injection, the 
synchrotron break frequency must be above the CHANDRA 
band, but this implies a magnetic field B < lOfiG. Even as- 
suming Tje — 1, with this small field all the models under- 
produce the radio emission by at least a factor 5. We find 
that it is possible to fit the synchrotron spectrum, including 
the average slope in the X-rays, only if 2.2 ^ 72 ^ 2.4. 

As to the other parameters, we notice the following: the 
value of rjM depends on age, with younger systems (i ~ 450 
yr) requiring a field B -^ 30/iG, and older ones (t ~ 650 



yr) requiring B ~ 20fiG, but it does not depend on 72; the 
value of He depends on 72, ranging from 4000 for 72 — 2.4 
to 4 X 10* for 72 = 2.2, but it is independent of age. All 
the models require a very high injection efficiency 77e ~ 1, 
to reproduce the INTEGRAL points. In Fig. |4l we show the 
spectrum derived according to the model in Tab. [T] 

The Lorentz factor corresponding to tc ranges between 
2x 10^ and 1.5x 10^, very close to the values derived in Crab. 
The one corresponding to the model of Fig. |3] is 7 x 10^. 
The minimum Lorentz factor is lower than ~ 5 x lO**, and 
the average Lorentz factor is ~ 7 x 10*, with an associated 
multiplicity that has to be greater than 10^. 

Kes 75 has been detected in gamma ra ys by HESS 
(JTerrier et al.l I2OO8I : iDiannati-Atai et all [JOOsI ') . The parti- 
cles responsible for TeV gamma-ray emission, in a purely 
leptonic model, are the ones emitting synchrotron radiation 
in the CHANDRA band. However the photon index in the 
TeV range is found to be ~ 2.3, far steeper than the average 
X-ray photon index of 1.9. Inverse Compton scattering on 
local synchrotron radiation gives a negligible contribution of 
TeV gamma-rays. Despite this being a young system, PSR 
J1846-0258 has a low spin-down luminosity, which results 
in a particle content about 2 orders of magnitude lower 
than in Crab and a correspondingly lower energy density 
of the synchrotron photon field. Inverse Compton scatter- 
ing on CMB photons is about 50 times stronger. However, 
even this contribution is about a factor of 10 below the ob- 
served fiux. In addition, scattering on the CMB is expected 
to occur in the Thompson regime and cannot account for 
the steepening of the spectral index in gamma-rays. A sug- 
gestion that has been made is that this steeper spectrum 
hints at a much warmer radiation background, for which IC 
should take place in the Klein-Nishina regime. This require- 
ment places the average black-body temperature of the seed 
background photons in the range 1000-2000 K. 

With a black-body spectrum at 1000 K, suppressed by 
a factor 5 x 10~®, it is possible to reproduce the correct 
gamma-ray luminosity. This contribution corresponds to a 
local enhancement of the infrared background which is only 
a factor of a few above the galactic average. Interestingly, as 
in the case of B1509, also Kes 75 is surrounded by a bright 
SNR shell, that might be responsible for the enhancement. 

The question arises if the TeV gamma-rays can be ex- 
plained by a hadronic component, and in particular by tt° 
decay. Given the very low magnetization, inferred from the 
synchrotron spectrum, it can be shown that only weak con- 
straints can be put on rjp-, in particular one can obtain a 
good fit of the radio and CHANDRA X-ray data, using 
rje ~ 0.5 (which however undepredicts INTEGRAL data) 
and a slightly higher nebular field, without violating the en- 
ergetics. However even by assuming half of the spin-down 
energy goes into protons, in order to have a significant con- 
tribution to the emission in the lOTeV range, at least IOMq 
of target thermal protons are needed. This is far in excess of 
the swept-up ejecta mass ~ O.IMq, even if consistent with 
the total ejecta mass (in this case one needs to assume that 
protons escape from the PWN and interact with the SNR). 
The main problem in this case is to properly reproduce the 
observed TeV spectrum: our model predicts a peak at 10^* 
Hz, instead of lO'^^ Hz as is observed, which is due to our 
choice of injection energy for the protons, which is tied to 7^ . 
Fermi observations should be able to distinguish between a 
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leptonic and a hadronic model, by constraining the emission 
below ITeV. To properly investigate whether n° decay is a 
viable possibility to explain the gamma-ray data, one would 
need a model for the diffusion of protons outside the nebula, 
which at the moment we do not have. The simplest possible 
estimate of the diffusion time of protons outside the nebula, 
and in the ejecta, gives values that are about an order of 
magnitude larger than the age of the nebula. 



6.5 Old Objects 

In this section we discuss two relatively old objects. The 
code has been developed for the investigation of systems also 
beyond the free expansion phase. However the late dynam- 
ics, especially if the pulsar kick velocity is important, can 
be quite complex. In addition old systems are usually very 
poorly constrained from the observational point of view. The 
age is generally known only as an order of magnitude esti- 
mate, and quite often either the central pulsar is not ob- 
served (as in the case of G327.1 and IC443), or the SNR is 
not observed. 

The following discussion highlights the difficulties that 
one has to face when trying to model sources for which both 
the quantity and the quality of the data are really poor. 



6.5.1 W44 

The SNR W44 is known to contain an old PWN, associated 
with pulsar PSR B1843-I-01. This is a 267 ms pulsar, with 
a characteristic dipole age of 2 0, 380 yr, and a spin -down 
luminosity of 4.3 x 10^^ erg s"^ l|Wolszczan et allll99ll '). The 
value of the braking index is not known, w hile the distance 
of th e pulsar is estimated to be 3.1 kpc (|Wolszczan et all 
[l992) for a typical electron density 0.03 cm~^. 

The SNR W44 (30 392) has an elongated shape with 
axes 25' x 11'. Its distance is estimated from HI absorp- 



tion to be 2.6 kpc, correspo nding to a typical SNR radius of 
11 — 13 pc (ICox et al.lll999l ). The small discrepancy between 
the estimated distance to the PSR and to the SNR is within 
the uncertainties, and the association is considered a secure 
one. This is a post-reverberation system, for which pres- 
sure balance between the PWN and the SNR is important. 
The central pressure is estimat ed from X-ray o bservations 
to be ~ 1 4 X 10^ djyne cm~^ dCox et al.lll999l) . Modeling 
the SNR (|Cox et al.lll999l : IShelton et ailll999i ) leads to the 
following estimates for the relevant dynamical parameters: 



Est 



lO''^ erg, 



Me 



5Mq, Po ~ 6 cm ; for the expan- 



sion velocity of the SNR , HI emission gives Vfg 



., „ _ 150 km 

s " (|Koo fc Heileslll995l ). A couple of remarks are here in 
order: first, all existing models assume the spin-down time 
of the pulsar as the age of the nebula; second, the quoted 
value of Vfs is inferred from a HI emitting ring structure 
which does not trace exactly the SNR, and es timates of the 
remn ant speed can be as high as 330 km s"'^ IjKoo fc Heiled 
Il995h . 

The PWN is observed in radio (JFrail et all Il996l : 
iGiacani et al.l 119971 ) with a typical luminosity of 200 mjy. 
It has a distinctive cometary shape with the pulsar located 
at the tip of a protruding finger of emission. The neb- 
ula is thought to be a transition object from the spherical 
shap e of young systems to later bow-shock like morpholo- 
gies (|van der Swaluw et al.ll2004h . Within this picture it is 
assumed that the PWN has already been crushed by the 
reverse shock, and that the pulsar has been displaced by its 
proper motion with respect to the core of the radio emis- 
sion. X-rays have been detected both with CHANDRA and 
XMM-Newton from t he head of the cometary nebula in th e 
vicinity of the pulsar l|Petre et al.ll2002l : iHarrus et aLll2006l ). 
and the fiux is measu red to be ^ 2.7 x 10~^^ erg cm"'^ s~ ^ 
in the 2-10 keV band (JHarrus et al.ll2006l : [Petre et al.ll2002l ). 
The radio emission is about twice as extended as the X- 
rays, suggesting effective cooling in the nebula, for which, 
however, there is no indication from the spatial behavior 
of the spectral index: this is found to be ^ 2.2 ± 0.3 with 
no appreciable variations with the distance from the pulsar 
(Harms et al. 2006). No IR or 7-ray emission is detected. At 
low radio frequencies the PWN is too weak to be detected 
against the SNR. 

The PWN is quite weak compared to the pulsar spin- 
down energy, and appears also to be quite small in size, even 
if a correct determination of the volume is problematic, given 
the shape. At the distance of the SNR, a typical value for 
the volume is ~ 0.5 pc'^, corresponding to a radius ^ 1.5 
ly. It is interesting to notice that models of the PWN have 
often assumed a much smaller age < 5000 yr IjPetre et al.l 
120021 ) ■ in order to fit the observed spectrum, which clearly 
contradicts what has been used to model the SNR. However 
such models are usually developed assuming that the spec- 
trum is Crab-like (in terms of the ratio between the radio 
and X-ray emission, and in terms of the location of the cool- 
ing break), but as we have shown through our modeling of 
young systems things can be different in different objects. 

The major proble m in model i ng th is system, as was 
already recognized bv IPetre et al.l (|2002l ). is the relatively 
low radio luminosity compared to the energy content of the 
nebula, and to the integrated pulsar spin-down power. This 
suggests that probably the age of the system is much lower. 
As a general rule, a younger age and a braking index less 



14 A^. Bucciantini, at al. 



W44 Luminos 



10-' 



lo-' 



10^ 



lo-' 



10" 



10^ 



10= 




10'= 10" 

Frequency (Hz) 
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Dashed line is the IC-CMB, dash-dotted line is the IC on the 
local IR background (same as for Crab nebula). Dotted line is 
the IC-SYN. 



than 3 correspond to less injected energy. In particular, we 
find that a braking index close to 2 and an age of about 
10, 000 — 15, 000 yr are required to reproduce the observed 
radio emission. We also require that the average magnetic 
field in the nebula must be well below equipartition. Indeed 
we find that, in order to avoid overproducing the radio emis- 
sion, r]M has to be < 0.002, corresponding to a nebular field 
below 10 /iG, even if in principle this refers to the average 
magnetic field, while the value in the head can be higher. 

This small value gains some support from the lack of 
spectral steepening with distance, which suggests that cool- 
ing is not important. Indeed the difference in size between 
the radio and X-ray nebula could be explained by the pres- 
ence of a much stronger magnetic field in the head of the 
nebula than in the body, resulting in the suppression of X- 
ray emission in the latter. In Table[5]we list the values of the 
parameters used to produce the curve in Fig. (5] Values of 
the SNR parameters are not dissimilar from what has bee n 
used in the literature l|Cox et al.lll999l : IShehon et al.lll999l ). 
The value we find for the central pressure is ~ 1.6 x 10®, in 
agreement with X-ray observations, but the expansion speed 
is found to be higher, Vfs ~ 300 km s^^, as a consequence 
of the younger assumed age. 

It is evident that the quantity and quality of existing 
data do not allow us to constrain the PWN model as much as 
it is possible for younger systems. Moreover, the shape is far 
from spherical, and geometrical effects connected with the 
existence of a bow shock cannot be modeled correctly within 
out approach. However, by repeating the same analysis we 
did in the previous section, we find that the parameter /le 
in the model is 8 x 10~^, which corresponds to a Lorentz 
factor at e^ equal to 2 x 10^, while the the average wind 
Lorentz factor is found to be 7^, <, 10*, corresponding to a 
multiplicity k >, 10*. 

One thing to notice is that the spectrum barely extends 
to X-ray frequencies. As pointed out above, this is due to 
the low value of the magnetization we adopt. If the magne- 
tization in the head is higher, as could be expected due to 



compression in the bow shock, this might enhance the X-ray 
emission. 

Another striking feature is the bump at around 10^* 
Hz: this is a typical feature of post-reverberation systems. 
During the reverberation phase, there is a typical energy 
at which synchrotron losses are balanced by adiabatic com- 
pression gains. Particles tend to accumulate at this energy, 
and this causes a bump in the spectrum, which survives at 
later times. One final point of interest is the fact that for 
this system (and for relatively old systems in general) IC 
scattering on the CMB gives a TeV flux that is comparable 
with that in the X-rays. 

6.5.2 K2/3 Kookaburra 

The "Kookaburra" is a complex of compact and extended 
radio/X-ray and 7-ray sources, that spans about one square 
degree along the galactic plane. A large circular thermal shell 
with a broad wing in the North-East and a narrower one in 
the South- West is revealed from radio images. Diffuse X-ray 
emission and point li ke sou rces have been detected by ASCA 
(Roberts et al. 19991 . l200ll ). XMM- Newton and CHANDRA 
i Ng et aLll2005^ . We are here interested in the North-East 
wing where at radio and X-ray frequencies a nebula is found, 
hosting a young and energetic pulsar, PSR J1420-6048. The 
pulsar location is also coincident with Te V 7-ray emission 
detected by HESS (| Aharonian et al.ll2006l ). 

PSR J1420-6048 is a 68.2 ms pulsar, with a character- 
istic dipole a ge of 13, 050 yr and a spin-down luminosity of 



10"" erg s~^ (JD'Amico et al.l l2001). The braking index is not 
known: as for the case of W44, this implies that the true age 
could differ substantially from the characteristic dipole age. 
The distance of the pu lsar is estimated to be 5.5±0.8 kpc 
(j Aharonian et al.ll200^ . 

Identifying the PWN associated with PSR J1420-6048 
is, however, rather problematic. The wing-like structure 
hosting the pulsar is usually referred to as K2. In coinci- 
dence with the PSR an enhancement of radio emission, usu- 
ally referred to as K3, is also observed. K2 has an extent of 
~ 15' x 10', a total flux at 20 cm of ~ 1 Jy, and a spectral 
index 0.2 ± 0.2, while K3 has an extent of about 3', a to- 
tal excess flux at 20 cm of ~ 2 mjy, and a spectral index 
0.4 ± 0.5 iRoberts et al.lll999l ). No IR detection has been 
reported. 

X-rays have been detected by ASCA in the 2-10 keV 
band, with an extent of about 7', an integrat ed flux of 4.8 x 
10~^ ^ erg cm~'^ s~^, and a spectral index 1.4 (jRoberts et al.l 
I2OOH ). A much more compact (0.5') but extended source has 
been detected by CHANDRA around the pulsar, with a total 
flux, extrapolated to the 2-10 keV band, of 1.3±0.1 4x 10"^^ 
erg c m~^ s~^, and with a spectral index 2.3 ± 0.9 (jNg et al.l 
I2OO5I ). 

Gamma-rays were flrst detected by EC RET 
]l99d), and more recently by HESS 



( Thompson et al 



(.Aharonian et al. 



[2OO6I ). The total HESS luminosity in the 
0.4-20 TeV band is 5.1 x 10^" erg s"\ for the assumed 
distance, and the spectral index is found to be ~ 2.2. The 
2a angular extension of the gamma-ray emission (~ 7') 
implies a nebular size ~ 11 pc. 

One immediately realizes that from a dynamical point 
of view these data are somewhat confusing. In particular we 
think that the standard interpretation, according to which 
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Table 5. Model parameters for W44 and K2/3 corresponding to the curves in figures [5] and |6] Units are the same as in Tab.[T] 
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the s mall K3 reRJon i s the PWN, while the K2 region is the 
SNR [Roberts et ahl (|200ll ). seems unlikely. Indeed, within 
the framework of our model, this assumption leads to un- 
physical values for the supernova energy and ejecta mass 
and to overpredict the radio emission. 

The larger extent of the ASCA source, compared to 
the CHANDRA one, goes against the expectation that the 
nebular size should shrink at higher frequencies because of 
synchrotron cooling. The difference in ASCA and CHAN- 
DRA spectral index might be due to uncertainties on the 
assumed value of Nh- The spectral index of the TeV emis- 
sion (interpreted as IC-CMB) is consistent with the steeper 
CHANDRA spectrum, but not with the flat ASCA one. 

Even the radio data are not conclusive: if one assumes 
that the PWN also contributes to the K2 emission, rather 
than just to K3, the radio emission changes by about 2 or- 
ders of magnitude. The uncertainties in the age, size and 
interpretation of the observed fluxes, in addition to the lack 
of a well defined SNR, seriously hamper the modeling of the 
system. Indeed, as for W44, we will present here just a sim- 
ple model fit to the data, without going into any detailed 
investigation of the parameter space. 

Again, as in the previous case, the major difficulty 
in accounting for what is observed relates to the low effi- 
ciency both in X-rays and radio, together with the ratio 
Lj/Lx > 1, which is suggestive of a large nebular content 
of low-energy particles. As in the case of W44 this implies 
that a decent fit to the data can only be achieved if one 
assumes a relatively young system, with a steep high en- 
ergy spectrum (similar to 3C58) and a weak magnetization 
(comparable with that in Kes75). Fig. [5] shows the result of 
the model corresponding to the parameters in Tab. \5\ The 
nebular radius is found to be ~ 5 pc in agreement with the 
HESS size. The magnetic field is inferred to be ^ 5/iG. It is 
clear that it is not possible to reproduce simultaneously the 
CHANDRA and ASCA data. To get the steeper CHANDRA 
spectrum one needs to assume a cooling break around lO'^® 
Hz. This is also necessary in order to repr oduce the correct 
TeV s pectrum. As already suggested by lAharonian et al.l 
l|2006l) , in the TeV range the IC-CMB dominates over the 
IC on the standard galactic background. Indeed our model 
under-predicts the TeV emission by a factor 4-5. The local 
photon background must be higher if the observed TeV flux 
has to be explained as IC emission: the TeV fiux shown in 
Fig. |6] was obtained by assuming an IR background in the 
form of Black-Body at 200K, with a photon density about 
a factor 4-5 higher than the standard galactic background. 
Given the complexity of the Kookaburra region, one cannot 
exclude such possibility. 

The parameter fie in the model is found to be 5 x 10""", 
which corresponds to a Lorentz factor at ec equal to 2 x 
10^, while the the average wind Lorentz factor is 7™ < 10*, 
corresponding to a multiplicity k > 10^. 




10'^ Mr 

Frequency (Hz) 

Figure 6. K3 /2 Kookaburra. Dat a points are: radio fro m from 
iRoberts et"all j figQaj . X-ray from [Roberts et al.l 1 I2OOII) . 7-ray 
from iNe et al.l (l2005t) : lAharonian et al.l ll2006h . Solid line is the 
total luminosity. Dashed line is the IC-CMB, dash-dotted line is 
the IC on the local IR background (assumed as a balck body at 
200K). Dotted line is the IC-SYN. 



7 DISCUSSION &: IMPLICATIONS 

Let us briefly review here our flndings, and what general 
conclusions can be drawn from the results of our attempt at 
modeling different objects, both young and old. 



7.1 Summary of the results 

First of all, it is remarkable that our one-zone model, de- 
spite its simplifications, has proven able to account for the 
observed Spectral Energy Distributions in all cases. This is 
already an interesting result, since the assumption of effi- 
cient mixing that is at its base is likely a poor approxima- 
tion at high energies (in the X-ray band), where the fiow 
pattern and the magnetic field structure in the inner nebula 
are expected to play a dominant role. 

The agreement between the model flt and the data is 
extremely good in the case of Crab, where the unknown pa- 
rameters are reduced to a minimum (essentially the SNR 
properties). Our model seems to work rather well for young 
systems in general, but the quality of the data and the lim- 
ited spectral coverage lead to larger uncertainties. In the 
case of strongly magnetized sources, like 3C58 and the neb- 
ula associated to PSR B1509, we find that in order to re- 
produce the data our model requires an energy input in the 
PWN which exceeds the pulsar release by 20-30%, judging 
from the current spin-down power. This conclusion does not 
seem to be affected by uncertainties in the distance to the 
sources: moving them closer to (further from) us would lower 
(increase) the estimated brightness, but at the same time it 
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Table 6. Summary of the inferred lower limits on multiplicity. 
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3C58 


MSH 15-52 


Kes75 


W44 


K2/K3 


10" 


5 X 10^ 


3 X 10^ 


10^ 


10^ 


105 



would make them smaller and younger (larger and older), 
leaving the discrepancy between the radiated and accumu- 
lated energy almost unaltered. We think it likely that this 
discrepancy is due to the simplifying assumptions of the one- 
zone model, which become progressively more important at 
large magnetizations. In partic ular, as also show n by recent 
multidimensional simulations (IVolpi et al.l 120081 ) . the radio 
and X-ray emitting particles might sample different mag- 
netic fields. At the same time it is worth keeping in mind 
that the estimate of the PSR spin-down power is based on 
the assumed canonical moment of inertia for neutron stars: 
I = 10*^ g cm^. A 20% discrepancy in two systems is com- 
patible with present uncertainties. 

As one can easily see from Tab.[l] given the already high 
value of the quantity ?7e -|- ??m, there is little energy (~ 20% 
at most) left to be carried by any other higher energy com- 
ponent (ions or leptons in the return current). The idea that 
a particle population with a larger Larmour radius could be 
carrying most of the wind energy was made attractive by 
the following three main reasons: these particles could help 
solve the problem of particle acceleration at the pulsar wind 
termination shock; they could explain the wisps variability 
in the Crab Nebula; finally, in the case of ions, they could 
solve the discrepancy between the predictions of 1-d MHD 
models and the gamma-ray flux observed from Crab. 

As far as particle acceleration is concerned, at a 
relativistic pair shock, efficient acceleration is shown to 
be p ossible only if th e shock magnetization is extremely 
low (ISpitkovskvlJ2008l) The 2-D MHD m odels of PWNe 
l|Del Zanna et al.l 120061 : ICamus et al.l 12003) . use a simple 
model of field reversal across the equator where a current 
sheet occurs, such that the appropriate conditions for accel- 
eration are thought to be realized only in a small sector of 
the pulsar wind termination shock, with no more than a few 
percent of the total wind energy flux flowing through it. 

On the other hand, the presence of an energetically sig- 
nificant ion component in the pulsar wind is found to lead 
to efficient acceleratio n of the pairs IjHoshino et al.l 1 19921 : 
lAmato fc Aronal200g) . This effect is the consequence of the 
large Larmor radius of the ions, whose gyration introduces 
long wavelength turbulence through cyclotron instability. 
The same result could come from the presence of high energy 
leptons, accelerated as a consequence of runaway dynamics 
in the equatorial current layer (Arons, in preparation). 

As to the wisp variability, Spi tkovskv fc AronsI (1200J) 
showed that this could be explained as the result of compres- 
sion waves associated with the high energy current carriers' 
(ions, in their case) gyration in the shock region: in order to 
reproduce the observed brightness contrast, again one would 
need the current carriers (the "beam") to be energetically 
significant: ^beammbeam / ^■wni± > 2k±. An alternate expla- 
nation of the wisp variability, which does without the kinetic 
effects of the high ener g y par ticles beam, has been recently 
shown bv lCamus et al.l 1 20091 ). who have proved that it can 



be recovered also within the framework of pure MHD, due 
to global instabilities of the termination shock. 

One final consideration concerns the possible detection 
of signatures of high energy protons. Modeling of the TeV 
emission from the Crab Nebula, based on the ID Kennel & 
Coroniti flow (jAtovan fc AharoniarJI 19961 ) was shown to un- 
derpredict the observed flux by a factor ~ 5, leading to sug- 
gest a possible contribution from the decay of neutral pions 
produced in nuclear collisions of relativistic protons. More 
recently, the presence of relativistic protons as the source of 
high- energy gamma- ray emiss ion was suggested also in the 
Vela PWN (|Horns et ffl/l2006f ). In the case of Vela, however, 
a later determination of the density of the thermal protons 
(that would serve as a target for nuclear collisions) resulted 
in too low a value and led to strongly questio n the initial 
claim (jLaMassa et aZ.1 120081 ). At the same time. IVolpi et al.l 
(J2008r ) showed that using a multidimensional model for the 
flow structure in Crab, the estimated IC-TeV emission is eas- 
ily overproduced, for the same parameters that allow one to 
better fit the lower energy synchrotron emission. This sug- 
gests that the results are indeed strongly model-dependent 
and that the discrepancies are likely due to problems with 
the adopted MHD model. 

This is also the conclusion we reach in this work. We can 
reproduce the TeV emission from Crab, and find results in 
agreement with Fermi observations without any need for a 
proton contribution. In the case of both B1509 and 3C58, our 
pair-|-magnetic field energetic is already exceeding by 20% 
the estimated PSR energy input, and fitting the gamma- 
rays through tt" decay requires that protons alone carry far 
more energy than the nebula is currently estimated to store. 
Finally, in the case of Kes75, in order to ascribe the excess 
TeV emission to protons, one needs to assume that protons 
can effectively diffuse out of the nebula and experience the 
whole ejecta mass as a target. 

Our model seems to constrain with good accuracy the 
value of £c- It is not possible to vary the value of tc by more 
than a factor of a few still satisfying the overall energet- 
ics and reproducing the radio and X-ray data. At the same 
time, our model allows us to derive, for each object, an up- 
per limit on the wind Lorentz factor and a lower limit on 
the pair multiplicity. The limits we put on 7^, and k are 
strictly valid only if the wind has a unique Lorentz factor 
and the energy scales em and tc respect the assumed scaling 
with time, proportional to the pulsar voltage. It is not easy 
to predict how a possible latitude dependence of the wind 
Lorentz factor, suc h as that included in ax isymmetric mod- 
els of PWNe {e.g. iDel Zanna et al.l l|2004l ')'). or a different 
dependence on time of the energy scales would affect our 
conclusions. 
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7.2 Interpretation of spectral breaks 

Our study has reinforced a long standing puzzle: the elec- 
tron (and positron) accelerator at work in PWNe knows 
how to create a spectrum convex in energy space, best 
re presented by a broken po wer law. As was advanced 
bv lKennel fc Coronitil (|l984br ). the spectral steepening ob- 
served between optical/soft X-ray energies and the harder 
(e > 10 keV) spectrum is well understood as the effect of 
synchrotron energy losses, with the pre-cooled spectrum a 
power law N{E) « E''''' , 2.1 < 72 < 2.8 from Table [T] 
The steepening between mid-infrared and optical, when ob- 
served, requires energy space structure in the accelerator 
(an "intrinsic break"). The high energy spectrum can be at- 
tributed to diffusive acceleration at the termination shock, 
at least qualitatively - Fermi acceleration in relativistic very 
low sigma shock s in the test particle limit yields a spectrum 
N{E) oc E-^'^ iKeshet fc WaxmanllJOOsI '). The radio spec- 
trum requires a much harder distribution of the lower energy 
particles, N(E) oc E'''^ 1.2 < 71 < 1.7. 

The standard model iKennel fc Coronitil (|l984al lbl). de- 
veloped to account for the radiation from Crab at near in- 
frared and shorter wavelengths, assigns the conversion of 
the pulsar wind energy in the particle spectra that emit 
the observed synchrotron radiation to diffusive Fermi accel- 
eration at the wind termination shock. Diffusive shock ac- 
celeration always shows particle spectra with a Maxwellian 
at low energy, plus a power law supra-thermal tail at 
high energy, with the temperature of the Maxwellian set 
by the shock jump conditions (except when the acceler- 
ation of the tail is very efficient). These properties are 
well exhibited by Particle-in-Cell simulations of relativis- 
tic shocks in unmagnetized e plasmas ( or, if magnetized , 
in upstream quasi-paralle l flow geometry l|Spitkovskvll2008l : 
ISironi fc Spitkovskvll2009l ')'l. Our models identify the transi- 
tion between the soft, high energy spectrum and the hard, 
low energy spectrum as being at the e nergy e^, typically 
lO'^ - 10^m±c^ (see Table [T]). In the (iKennel fc Coronitil 
Il984bl ) model, this energy was identified as the lower cut- 
off of the shock accelerated power law, and that energy was 
identified as the "temperature" ~ 7u,m±c^ of the fiow down- 
stream of the termination shock, thus giving rise to the be- 
lief that pulsar winds have upstream fiow Lorentz factor 
7» ~ 10^ 

The Kennel and Coroniti model and its descendants de- 
liberately neglected the radio emission from the Crab Neb- 
ula (and, by extension, other PWNe.) The large (by number) 
population of particles with E < ec suggests 7^, to be much 
smaller, if the radio and mid- to far-IR emitting particles 
come from the pulsar. That the pulsars are the most likely 
source of the low energy particles in each of the nebulae gains 
support fr om the existence of the radio "wisp" features in 
the Crab (|Bietenholal2004l ) closely associated with the sim- 
ilar time variable structures seen in optical and X-ray imag- 
ing, as well as similar structure seen in 3C58 (|Bietenhola 
120061 ) - those structures are clearly coincident with the ter- 
mination shock. Then, in the regions of the fiow populated 
by the large particle fiux feeding the low energy population, 
7™ ~ 10* or smaller. 

The low energy particle spectrum definitely is not 
a single temperature relativistic Maxwellian, with T ~ 
Ec. Since the termination shock is not spherical, becom- 



ing more oblique in higher r otational latitude re gions 
(JKomissarov fc Lvubarskvll2004l : rDel Zanna et al.ll2004l '). the 
post-shock temperature declines with increasing latitude, 
suggesting the low energy particles enter in the polar regions 
of the outfiow. However, replicating the observed spectrum 
as the envelope of a sequence of Maxwellians requires a fac- 
tor of ~ 100 variation in temperature, which requires most 
of the mass fiux being nearest the rotation poles. There is 
no sign of pole to equator asymmetry in the radio emission 
near the pulsars, except for the radio wisp features in the 
Crab, which are components of the immediate post-shock 
fiow, as judged from the optical and harder photon emis- 
sion. However we cannot rule out that mixing and effects 
related to integration along the line of sight might prevent 
the detection of the implied variations. 

However, let us assume in the following, that also the 
low-energy particle spectrum results from some acceleration 
process, rather than from the convolution of different ther- 
mal distributions, and let us speculate on the nature of such 
a process. 

Turbulence and associated Fermi II acceleration in and 
around the termination shock is one appealing possibility. 
The "wisp" motio n, observed from the radio through X-rays 
(see|HesteiJ (|2008l ) for a comprehensive review) has recently 
been interpreted as the result of the strongly variable termi- 
nation shock structure foun d in high re solution MHD simu- 
lations of the Crab Nebula (|Camus et al. 2009 ). The shock 
instability implied is a termination shock variant of the 
Standing Accretion Shock Instability, with outer scale vari- 
able velocity 5v ^ v ^ Q.Qc and length scale (5r~r~0.5 — 1 
light years. The magnetized motions observed in the simula- 
tions (which do well in replicating the time variable spatial 
structure observed in the nebular "wisps ") can act as an 
accelerator through t he Fe rmi II process (|Kardashevlll962l : 
[Stawarz & PctrosianI I2OO8I ) . creating electron and positron 
spectra N{E) = NoiE/Eo)-'. 

In quasi-linear models of Fermi II acceleration 
in isotropic small amplitude Alfven turbulence {e.g. 
IStawarz fc PetrosianI (|2008[ )). s = y^9/4 -I- e--l/2, where e = 
Taccei/Teacape. Thcsc modcls iuclude an analogue of scatter- 
ing from large amplitude magnetic inhomogeneities (Fermi's 
original model), in the case where the wave energy spectrum 
scales in proportion to k~^ , with k the wavenumber. This 
case is germane to acceleration due to pair interaction with 
large scale moving magnetic fiuctuations: the acceleration 
time scale for interaction with large scale moving "eddies" 
is Taccei ~ Xo/c{c/veddy)^iB/5B)^ ~, whcrc Ao is the outer 
scale of t he turbulence, com parable to the shock radius in 
our case (jCamus et al.ll2009l ) . The high downstream fiow ve- 
locity suggests particles escape the turbulence zone through 
advective loss rather than diffusive escape (at energies much 
less than a PeV, microscopic diffusion across B, even at the 
Bohm rate, is negligible compared to the advective losses). 
Then at fiow speed ~ c/3, Tescape ~ SRshock/c, leading to 
e ~ {l/3){B/SB)'^{Xo/Rshock)- Taking the parameters to be 
unity yields e '^ 1/3, in which case the accelerated spectrum 
is very hard: N{E) oc E~^'^ which is similar to the spectrum 
inferred for the radio emitting pairs in PWNe. Because of 
the strong radial mixing observed in the MHD models of 
PWNe, the spectrum created in this turbulent acceleration 
zone will fill the whole body of the synchrotron emitting neb- 
ulae, which provides a natural explanation of the observed 



A^. Bucciantini, at al. 



lack of gradients in the radio spectral index of the Crab 
l|Bietenholall992[). Thus, we rev i ve, in a modern form, the 
suggestion of iBarnes fc Scargld ()l973l ) , that the wisp mo- 
tions are responsible for a part of the particle acceleration 
required to account for PWNe synchrotron emission. 

Turning this idea into a full physical model may re- 
quire extending the MHD models to three dimensions. All 
models to date have been 2D, axisymmetric, with exclu- 
sively toroidal magnetic field. Particles then cannot scatter 
radially, in order to sample adjacent moving eddies, but are 
confined to toroidal flux tubes, which move in and out in 
radius around average positions in the outflow (and inflow, 
at higher latitudes) in the turbulent region. It is possible 
that magnetic pumping in these flux tubes might substitute 
for particle scattering in qua si-isotropic tur b ulence as an ac- 
celera tion mechanism (e.g.. iMelrosd (| 19691 ) : iKuiipers et al\ 
119971)). However, polarization studies already demonstrate 
that there are substantial poloidal magnetic fields in PWNe 
inner regions, at least i n the Crab. The ter mination shock in- 
stability identified by IjCamus et al.ll2009l ) is not likely to be 
restricted to 2D - perhaps the assumption of quasi-isotropic 
turbulence, with forcing scale ~ Rshock, is the most natural 
starting point for further pursuit of this idea. 

It is clear that a key point for the model to succeed is 
that the proposed turbulently accelerated spectrum at low 
energy merges smoothly into the shock accelerated spectrum 
at high energy. Since the macroscopic turbulence is closely 
associated with the shock, a smooth merger is at least think- 
able. In this context the meaning of ec changes completely 
with respect to the Kennel and Coroniti interpretation: this 
should not be identified with the post-shock temperature 
(which according to the present model is about an order of 
magnitude lower), but rather with the energy above which 
the low energy accelerator (Fermi II in the turbulent flow, 
in the suggestion made here) fades out and relativistic DSA 
takes over. 



7.3 How Do We Move Forward? 

Multi-D high resolution MHD and PIC simulations of the 
termination shock region separating the freely expanding 
wind from the subsonically expanding nebulae will shed light 
on the possibly turbulent flow, and on whether that turbu- 
lence can act as the low energy accelerator discussed previ- 
ously. A useful first step would be to couple the diffusion- 
advection equation in energy space to the MHD calculations. 
Observationally, infrared observati ons of t he young PWNe, 
using the Herschel space telescope (jPilbrat t 2009) as well as 
near infrared instruments on the ground and in space, will 
greatly improve our understanding of the physics behind 
the broken power laws in energy space inferred here for the 
injected particle distributions. These instruments will have 
more than sufficient angular resolution to resolve the PWNe, 
thus emphasizing the need for multidimensional models to 
quantitatively interpret the observations. 

For all the objects studied, the plasma injection rates, 
here inferred from the evolutionary models, not simply taken 
from the average over the nebular lifetimes, exceed 10^ times 
the Goldreich- Julian rate. Such high rates can be under- 
stood only as the result of pair creation within the pulsars' 
magnetospheres - given the spectral continuity in the SEDs, 
models which rely on the radio and far-IR emitting parti- 



cles, which dominate the injection rates, being accelerated 
from the thermal plasma in and around the nebulae are less 
likely than injection from the pulsars . Our resu lts support 
and extend a long held suspicion (Gallant et 0^.112002 ) . pre- 
viously based only on analysis of the Crab Nebula, that the 
particle loss rates from young, high voltage pulsars are sub- 
stantially in excess of the inferences of particle outflow from 
all known magnetospheric pair creation models. The answer 
may lie in intrinsic ti me dependence of pair creation within 
pulsar polar caps fe.o. iTimokhinI (|2009|)), or in the so far un- 
explored outer magnetosphere accelerators associated with 
the boundary la yer return cu rrents separating the closed and 
open field lines ( Arong[2009|) . Fermi gamma ray observations 
are and will be useful in constraining magnetospheric pair 
creation models. 
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